🌳 Roadside Deforestation & Resilience: A Panarchy Framework¶
Assessing Socio-Ecological Impacts in Charghat using PRA & Geospatial Python¶
📝 Overview¶
This project investigates the severe socio-ecological crisis triggered by the widening of the Z-6006 Highway in Charghat, Rajshahi, which led to the removal of over 4,000 mature roadside trees.
Using the Panarchy Adaptive Cycle as a theoretical lens, the study integrates Participatory Rural Appraisal (PRA) with Geospatial Python (Google Earth Engine) to quantify the transition from a stable "Conservation" phase to a chaotic "Release" phase. The project assesses the cascading impacts—heat stress, air pollution, and livelihood loss—and co-designs a Locally Led Adaptation (LLA) framework to guide the system toward resilient "Renewal".
🎯 Objectives¶
- Objective 1: Assess physical, socio-economic, and environmental impacts using the Panarchy Theory (Conservation $\to$ Release $\to$ Reorganization).
- Objective 2: Validate community perceptions of heat and pollution using Remote Sensing (LST, NDVI, AOD).
- Objective 3: Develop a community-driven Nature-Based Solutions (NBS) framework for future resilience.
🛠️ Tools & Technologies¶
-lightgrey)
⚙️ Methodology¶
| Step | Description |
|---|---|
| 1. Data Collection | Primary: PRA tools (Social Mapping, FGDs, Problem Trees). Secondary: Satellite data via GEE (Sentinel-2, Landsat 8, MODIS, CHIRPS). |
| 2. Preprocessing | Used Python (geemap) to filter cloud cover, harmonize Sentinel-2 data, and compute indices (NDVI, AOD) clipped to the Charghat ROI. |
| 3. Analysis | Qualitative: Panarchy phase interpretation of community narratives. Spatial: Time-series analysis of LST (2016-2024) and Vegetation loss. |
| 4. Visualization | Created Sankey Diagrams to map causal flows of the crisis and Heatmaps to quantify stakeholder priorities. |
| 5. Validation | Triangulated PRA insights (e.g., "it feels hotter") with quantitative LST data (e.g., +5°C increase). |
Presentation¶
Total Pages: 31 --- Page 1 ---
--- Page 2 ---
--- Page 3 ---
--- Page 4 ---
--- Page 5 ---
--- Page 6 ---
--- Page 7 ---
--- Page 8 ---
--- Page 9 ---
--- Page 10 ---
--- Page 11 ---
--- Page 12 ---
--- Page 13 ---
--- Page 14 ---
--- Page 15 ---
--- Page 16 ---
--- Page 17 ---
--- Page 18 ---
--- Page 19 ---
--- Page 20 ---
--- Page 21 ---
--- Page 22 ---
--- Page 23 ---
--- Page 24 ---
--- Page 25 ---
--- Page 26 ---
--- Page 27 ---
--- Page 28 ---
--- Page 29 ---
--- Page 30 ---
--- Page 31 ---
Codes ⬇️¶
Libraries, Authentication & Study Area¶
import ee
import geemap
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Initialize GEE
try:
ee.Initialize()
except Exception:
ee.Authenticate()
ee.Initialize()
# Load Study Area
roi = ee.FeatureCollection('projects/geo-gee-imtiaj/assets/Charghat_upazila_Rajshahi')
print("Setup Complete. ROI Loaded.")
Setup Complete. ROI Loaded.
VEGETATION ANALYSIS (NDVI - Sentinel 2 10m)¶
# --- 1. NDVI ANALYSIS & EXPORTS ---
print("Processing NDVI...")
# Function to get NDVI
def get_ndvi(year):
s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED') \
.filterBounds(roi).filterDate(f'{year}-03-01', f'{year}-05-31') \
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)).median().clip(roi)
return s2.normalizedDifference(['B8', 'B4']).rename(f'NDVI_{year}')
# A. Generate Images
ndvi_16 = get_ndvi(2016)
ndvi_20 = get_ndvi(2020)
ndvi_24 = get_ndvi(2024)
ndvi_change = ndvi_24.subtract(ndvi_16).rename('NDVI_Change')
# B. Calculate Stats (Full Time Series 2016-2024)
years = list(range(2016, 2025))
stats_data = []
for y in years:
img = get_ndvi(y)
# 1. Mean NDVI
mean = img.reduceRegion(ee.Reducer.mean(), roi, 30, maxPixels=1e9).get(f'NDVI_{y}').getInfo()
# 2. Area Calculation (Vegetation > 0.3)
# Create a mask where NDVI > 0.3 (Vegetation)
veg_mask = img.gt(0.3)
# Calculate area of masked pixels
area_img = veg_mask.multiply(ee.Image.pixelArea())
area_sqm = area_img.reduceRegion(ee.Reducer.sum(), roi, 30, maxPixels=1e9).get(f'NDVI_{y}').getInfo()
area_ha = area_sqm / 10000 if area_sqm else 0 # Convert m2 to Hectares
stats_data.append({'Year': y, 'Mean_NDVI': mean, 'Veg_Area_Ha': area_ha})
df_ndvi = pd.DataFrame(stats_data)
# C. Visualization (Charts)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
# Line Chart (Start to End Change)
sns.lineplot(data=df_ndvi, x='Year', y='Mean_NDVI', ax=ax1, color='green', marker='o')
ax1.set_title('NDVI Trend (2016-2024)')
ax1.set_ylabel('Average NDVI')
ax1.grid(True)
# Bar Chart (Area Change)
sns.barplot(data=df_ndvi[df_ndvi['Year'].isin([2016, 2020, 2024])], x='Year', y='Veg_Area_Ha', ax=ax2, palette='Greens')
ax2.set_title('Vegetation Area Change (Hectares)')
ax2.set_ylabel('Area (ha)')
plt.tight_layout()
plt.savefig('NDVI_Analysis_Charts.png')
plt.show()
# D. Map Display (Side by Side)
Map1 = geemap.Map()
Map1.centerObject(roi, 12)
vis = {'min': 0, 'max': 0.7, 'palette': ['red', 'yellow', 'green']}
change_vis = {'min': -0.3, 'max': 0.3, 'palette': ['red', 'white', 'green']}
Map1.addLayer(ndvi_16, vis, 'NDVI 2016')
Map1.addLayer(ndvi_20, vis, 'NDVI 2020')
Map1.addLayer(ndvi_24, vis, 'NDVI 2024')
Map1.addLayer(ndvi_change, change_vis, 'NDVI Change (16-24)')
display(Map1)
# E. Export to Drive
geemap.ee_export_image_to_drive(ndvi_16, description='NDVI_2016', folder='Charghat_GIS', scale=10, region=roi.geometry())
geemap.ee_export_image_to_drive(ndvi_24, description='NDVI_2024', folder='Charghat_GIS', scale=10, region=roi.geometry())
geemap.ee_export_image_to_drive(ndvi_change, description='NDVI_Change', folder='Charghat_GIS', scale=10, region=roi.geometry())
print("NDVI Images exported to Drive folder 'Charghat_GIS'.")
Processing NDVI...
Map(center=[24.30210585383211, 88.77349383046975], controls=(WidgetControl(options=['position', 'transparent_b…
NDVI Images exported to Drive folder 'Charghat_GIS'.
TEMPERATURE ANALYSIS (LST - Landsat 8 30m)¶
# --- CELL: LST ANALYSIS (Temperature) - FIXED ---
import ee
import geemap
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Initialize (if not already done)
try:
ee.Initialize()
except Exception:
ee.Authenticate()
ee.Initialize()
# 1. Setup
roi = ee.FeatureCollection('projects/geo-gee-imtiaj/assets/Charghat_upazila_Rajshahi')
years = list(range(2016, 2025)) # Define years explicitly
print("Processing LST (Temperature)...")
def get_lst(year):
# Filter: March-May
# FIX: Increased CLOUD_COVER to 60 to ensure we find images in cloudy years
l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') \
.filterBounds(roi) \
.filterDate(f'{year}-03-01', f'{year}-05-31') \
.filter(ee.Filter.lt('CLOUD_COVER', 60)) \
.median() \
.clip(roi)
# Calculate LST (Kelvin to Celsius)
# Band 10 is Thermal Infrared Sensor (TIRS)
lst = l8.select('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15).rename(f'LST_{year}')
return lst
# A. Generate Images for Display
# Note: If 2016 is too cloudy, the map layer might be gray/empty.
# In that case, try changing the year to 2017.
lst_16 = get_lst(2016)
lst_24 = get_lst(2024)
lst_change = lst_24.subtract(lst_16).rename('LST_Change')
# B. Calculate Stats (Full Time Series) with Error Handling
lst_stats = []
print("Extracting statistics (this handles cloudy years)...")
for y in years:
try:
img = get_lst(y)
# Attempt to get the mean
mean = img.reduceRegion(ee.Reducer.mean(), roi, 30, maxPixels=1e9).get(f'LST_{y}').getInfo()
# Check if mean is None (empty image)
if mean is not None:
lst_stats.append({'Year': y, 'Mean_LST': mean})
else:
print(f"Skipping {y}: No clear images found.")
except Exception as e:
print(f"Error processing {y}: {e}")
df_lst = pd.DataFrame(lst_stats)
# C. Visualization (Charts)
if not df_lst.empty:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
# Line Chart (Trend)
sns.lineplot(data=df_lst, x='Year', y='Mean_LST', ax=ax1, color='red', marker='s')
ax1.set_title('Summer Temperature Trend (2016-2024)')
ax1.set_ylabel('Temp (°C)')
ax1.grid(True)
# Bar Chart (Comparison)
# Only plot years that actually exist in the dataframe
available_years = [y for y in [2016, 2020, 2024] if y in df_lst['Year'].values]
sns.barplot(data=df_lst[df_lst['Year'].isin(available_years)], x='Year', y='Mean_LST', ax=ax2, palette='Reds')
ax2.set_title('Avg Summer Temperature Comparison')
ax2.set_ylim(20, 40)
plt.tight_layout()
plt.savefig('LST_Analysis_Charts.png')
plt.show()
else:
print("No LST data could be extracted. Try widening the date range or cloud filter.")
# D. Map Display
# This visualizes the Urban Heat Island Effect concepts
Map2 = geemap.Map()
Map2.centerObject(roi, 12)
vis_lst = {'min': 20, 'max': 40, 'palette': ['blue', 'cyan', 'yellow', 'red']}
vis_diff = {'min': 0, 'max': 5, 'palette': ['white', 'yellow', 'red']}
Map2.addLayer(lst_16, vis_lst, 'LST 2016')
Map2.addLayer(lst_24, vis_lst, 'LST 2024')
Map2.addLayer(lst_change, vis_diff, 'LST Change (16-24)')
display(Map2)
# E. Export
print("Starting Exports...")
geemap.ee_export_image_to_drive(lst_16, description='LST_2016', folder='Charghat_GIS', scale=30, region=roi.geometry())
geemap.ee_export_image_to_drive(lst_24, description='LST_2024', folder='Charghat_GIS', scale=30, region=roi.geometry())
geemap.ee_export_image_to_drive(lst_change, description='LST_Change', folder='Charghat_GIS', scale=30, region=roi.geometry())
print("LST Export tasks sent to GEE.")
Processing LST (Temperature)... Extracting statistics (this handles cloudy years)...
Map(center=[24.30210585383211, 88.77349383046975], controls=(WidgetControl(options=['position', 'transparent_b…
Starting Exports... LST Export tasks sent to GEE.
Precipitation Analysis¶
# --- 4. PRECIPITATION ---
print("Processing Precipitation...")
rain_stats = []
for y in years:
img = ee.ImageCollection('UCSB-CHG/CHIRPS/PENTAD') \
.filterDate(f'{y}-01-01', f'{y}-12-31').sum().clip(roi)
val = img.reduceRegion(ee.Reducer.mean(), roi, 5000).get('precipitation').getInfo()
rain_stats.append({'Year': y, 'Rainfall_mm': val})
df_rain = pd.DataFrame(rain_stats)
# Chart
plt.figure(figsize=(10, 5))
plt.bar(df_rain['Year'], df_rain['Rainfall_mm'], color='#4a90e2')
plt.plot(df_rain['Year'], df_rain['Rainfall_mm'], color='navy', marker='o') # Line overlay
plt.title('Annual Total Precipitation')
plt.ylabel('Rainfall (mm)')
plt.savefig('Precipitation_Chart.png')
plt.show()
print("All tasks complete. Check Google Drive folder 'Charghat_GIS' for maps.")
Processing Precipitation...
All tasks complete. Check Google Drive folder 'Charghat_GIS' for maps.
Air Quality (NO2) Analysis¶
# --- 3. AIR QUALITY (Corrected using MODIS AOD) ---
import ee
import geemap
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Initialize
try:
ee.Initialize()
except Exception:
ee.Authenticate()
ee.Initialize()
# Study Area
roi = ee.FeatureCollection('projects/geo-gee-imtiaj/assets/Charghat_upazila_Rajshahi')
years = list(range(2016, 2025))
print("Processing Air Quality (Aerosol Optical Depth)...")
# Data Container
air_stats = []
for y in years:
try:
# DATASET: MODIS Terra Monthly Level 3 (Reliable 2000-Present)
# ID: MODIS/061/MOD08_M3
img = ee.ImageCollection('MODIS/061/MOD08_M3') \
.filterDate(f'{y}-01-01', f'{y}-12-31') \
.mean() \
.clip(roi)
# BAND: 'Aerosol_Optical_Depth_Land_Ocean_Mean_Mean'
# This specific band was confirmed from your previous error logs
aod_val = img.select('Aerosol_Optical_Depth_Land_Ocean_Mean_Mean') \
.reduceRegion(ee.Reducer.mean(), roi, 5000) \
.get('Aerosol_Optical_Depth_Land_Ocean_Mean_Mean') \
.getInfo()
if aod_val is not None:
air_stats.append({'Year': y, 'AOD': aod_val})
else:
print(f"No data for {y}")
except Exception as e:
print(f"Error in {y}: {e}")
# Visualization
if air_stats:
df_air = pd.DataFrame(air_stats)
# 1. Line Chart
plt.figure(figsize=(10, 5))
sns.lineplot(data=df_air, x='Year', y='AOD', marker='^', color='purple', linewidth=2.5)
plt.title('Air Pollution Trend (Aerosol Optical Depth)')
plt.ylabel('Mean AOD (Higher = Hazier/Dustier)')
plt.grid(True, linestyle='--')
plt.savefig('Air_Quality_Chart.png')
plt.show()
# 2. Export 2024 Map
# We generate a map for the latest year to show spatial distribution
print("Exporting 2024 AOD Map...")
air_img_24 = ee.ImageCollection('MODIS/061/MOD08_M3') \
.filterDate('2024-01-01', '2024-12-31').mean() \
.select('Aerosol_Optical_Depth_Land_Ocean_Mean_Mean') \
.clip(roi)
geemap.ee_export_image_to_drive(
air_img_24,
description='AOD_AirQuality_2024',
folder='Charghat_GIS',
scale=1000,
region=roi.geometry()
)
print("Air Quality Export sent to Drive.")
else:
print("No Air Quality statistics could be generated.")
Processing Air Quality (Aerosol Optical Depth)...
Exporting 2024 AOD Map... Air Quality Export sent to Drive.
import ee
import geemap
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import seaborn as sns
import numpy as np
import os
# --- 1. INITIALIZATION ---
try:
ee.Initialize()
except Exception:
ee.Authenticate()
ee.Initialize()
# Load Study Area
roi = ee.FeatureCollection('projects/geo-gee-imtiaj/assets/Charghat_upazila_Rajshahi')
print("Setup Complete. ROI Loaded.")
# Create a local output folder
out_dir = os.path.join(os.getcwd(), 'Charghat_Outputs')
if not os.path.exists(out_dir):
os.makedirs(out_dir)
print(f"All files will be saved to: {out_dir}")
# ==========================================
# HELPER FUNCTION: CREATE MAP LAYOUT (JPEG)
# ==========================================
def export_map_layout(images, titles, filename, cmap, min_val, max_val, label):
"""
Exports a professional side-by-side map layout to JPEG.
"""
print(f"Generating Map Layout: {filename}...")
# Fetch data as numpy arrays (Scale 100m for visualization speed)
# Note: High-res analysis is done in GEE; this is just for the JPEG map visual
data_list = []
for img in images:
try:
# Handle masking for clean borders
arr = geemap.ee_to_numpy(img, region=roi.geometry(), scale=100)
data_list.append(arr)
except Exception as e:
print(f"Error fetching map data: {e}")
return
# Create Plot
num_maps = len(images)
fig, axes = plt.subplots(1, num_maps, figsize=(6 * num_maps, 6), constrained_layout=True)
if num_maps == 1: axes = [axes]
# Plot each map
for ax, arr, title in zip(axes, data_list, titles):
# Handle NoData (which comes as extremely low/high numbers in numpy)
arr = np.nan_to_num(arr, nan=np.nan)
im = ax.imshow(arr, cmap=cmap, vmin=min_val, vmax=max_val)
ax.set_title(title, fontsize=14, fontweight='bold', pad=10)
# Neatline (Border)
for spine in ax.spines.values():
spine.set_edgecolor('black')
spine.set_linewidth(2)
# Remove ticks for clean map look
ax.set_xticks([])
ax.set_yticks([])
# North Arrow (Simple text representation)
ax.text(0.95, 0.95, 'N \u2191', transform=ax.transAxes,
ha='right', va='top', fontsize=16, fontweight='bold', color='black')
# Add Colorbar (Legend)
cbar = fig.colorbar(im, ax=axes, orientation='horizontal', fraction=0.05, pad=0.05)
cbar.set_label(label, fontsize=12)
# Save locally
save_path = os.path.join(out_dir, filename)
plt.savefig(save_path, dpi=300, bbox_inches='tight')
plt.close()
print(f"Saved Map: {save_path}")
# ==========================================
# 2. VECTOR EXPORT (Shapefile)
# ==========================================
print("Exporting Boundary Shapefile...")
geemap.ee_export_vector(roi, filename=os.path.join(out_dir, 'Charghat_Boundary.shp'))
# ==========================================
# 3. NDVI ANALYSIS & MAPS
# ==========================================
print("\n--- PROCESSING NDVI ---")
def get_ndvi(year):
s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED') \
.filterBounds(roi).filterDate(f'{year}-03-01', f'{year}-05-31') \
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)).median().clip(roi)
return s2.normalizedDifference(['B8', 'B4']).rename('NDVI')
# A. Generate Images
# ndvi_16 = get_ndvi(2016)
ndvi_20 = get_ndvi(2020)
ndvi_24 = get_ndvi(2024)
ndvi_change = ndvi_24.subtract(ndvi_20).rename('NDVI_Change')
# B. Export Raster Data (TIF)
geemap.ee_export_image(ndvi_20, filename=os.path.join(out_dir, 'NDVI_2016.tif'), scale=10, region=roi.geometry())
geemap.ee_export_image(ndvi_24, filename=os.path.join(out_dir, 'NDVI_2024.tif'), scale=10, region=roi.geometry())
geemap.ee_export_image(ndvi_change, filename=os.path.join(out_dir, 'NDVI_Change.tif'), scale=10, region=roi.geometry())
# C. Generate JPEG Maps (Side-by-Side)
export_map_layout(
images=[ndvi_20, ndvi_24],
titles=['NDVI 2020', 'NDVI 2024'],
filename='Map_NDVI_Comparison.jpg',
cmap='RdYlGn', min_val=0, max_val=0.6, label='Vegetation Index (NDVI)'
)
# D. Generate Change Map
export_map_layout(
images=[ndvi_change],
titles=['NDVI Change (2020-2024)'],
filename='Map_NDVI_Change.jpg',
cmap='RdYlGn', min_val=-0.3, max_val=0.3, label='Change (Red=Loss, Green=Gain)'
)
# ==========================================
# 4. LST ANALYSIS & MAPS
# ==========================================
print("\n--- PROCESSING LST ---")
def get_lst(year):
l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') \
.filterBounds(roi).filterDate(f'{year}-03-01', f'{year}-05-31') \
.filter(ee.Filter.lt('CLOUD_COVER', 60)).median().clip(roi) # Higher tolerance for clouds
return l8.select('ST_B10').multiply(0.00341802).add(149.0).subtract(273.15)
# A. Generate Images
# lst_16 = get_lst(2016)
lst_20 = get_lst(2020)
lst_24 = get_lst(2024)
lst_change = lst_24.subtract(lst_16)
# B. Export Raster Data (TIF)
geemap.ee_export_image(lst_20, filename=os.path.join(out_dir, 'LST_2016.tif'), scale=30, region=roi.geometry())
geemap.ee_export_image(lst_24, filename=os.path.join(out_dir, 'LST_2024.tif'), scale=30, region=roi.geometry())
# C. Generate JPEG Maps
export_map_layout(
images=[lst_20, lst_24],
titles=['LST 2020', 'LST 2024'],
filename='Map_LST_Comparison.jpg',
cmap='RdYlBu_r', min_val=22, max_val=38, label='Temperature (°C)'
)
export_map_layout(
images=[lst_change],
titles=['LST Increase (2020-2024)'],
filename='Map_LST_Change.jpg',
cmap='inferno', min_val=0, max_val=5, label='Temp Increase (°C)'
)
# ==========================================
# 5. PRECIPITATION & AIR QUALITY MAPS
# ==========================================
print("\n--- PROCESSING PRECIP & AIR ---")
# Precipitation Map (2024 Total)
precip_img = ee.ImageCollection('UCSB-CHG/CHIRPS/PENTAD') \
.filterDate('2024-01-01', '2024-12-31').sum().clip(roi)
geemap.ee_export_image(precip_img, filename=os.path.join(out_dir, 'Precipitation_2024.tif'), scale=1000, region=roi.geometry())
export_map_layout(
images=[precip_img],
titles=['Total Rainfall 2024'],
filename='Map_Precipitation_2024.jpg',
cmap='Blues', min_val=1000, max_val=2000, label='Rainfall (mm)'
)
# Air Quality Map (MODIS AOD 2024)
# Note: Using MODIS because OMI NO2 ID was broken in your region
air_img = ee.ImageCollection('MODIS/061/MOD08_M3') \
.filterDate('2024-01-01', '2024-12-31').mean() \
.select('Aerosol_Optical_Depth_Land_Ocean_Mean_Mean').clip(roi)
geemap.ee_export_image(air_img, filename=os.path.join(out_dir, 'AirQuality_AOD_2024.tif'), scale=1000, region=roi.geometry())
export_map_layout(
images=[air_img],
titles=['Air Quality (Aerosol Optical Depth) 2024'],
filename='Map_AirQuality_2024.jpg',
cmap='YlOrBr', min_val=0.3, max_val=1.0, label='AOD (Unitless)'
)
print(f"\nAll Tasks Completed Successfully.")
print(f"Check the folder: {out_dir}")
📊 Results & Insights¶
- 🌍 System Collapse: The removal of trees triggered a "Release" phase, causing a sharp decline in NDVI and a spike in Aerosol Optical Depth (AOD).
- 📈 Heat Island Effect: Land Surface Temperature (LST) analysis revealed that summer temperatures in 2024 reached up to 39°C, a significant increase from the 2016 baseline.
- 🧭 Livelihood Shock: Code co-occurrence analysis showed that Livelihood Insecurity was the highest priority (100% intensity) for farmers and market vendors, overshadowing biodiversity concerns.
🗂️ Data Sources¶
| Source | Description | Link |
|---|---|---|
| Sentinel-2 | Vegetation Index (NDVI) 10m | Copernicus |
| Landsat 8 | Land Surface Temperature (LST) 30m | USGS EarthExplorer |
| MODIS Terra | Aerosol Optical Depth (AOD) | NASA Earthdata |
| Primary Data | FGD, KII, and Social Mapping | Field Survey (Charghat, 2024) |
🔖 Tags¶
Participatory Planning Geospatial Python Panarchy Theory Google Earth Engine Climate Resilience LST Analysis Nature-Based Solutions Rajshahi